Systematic Parameterization of Ion–Surfactant Interactions in Dissipative Particle Dynamics Using Setschenow Coefficients

Dissipative particle dynamics (DPD) simulations of nonionic surfactants with an added salt show that the Setschenow relationship is reproduced; that is, the critical micelle concentration is log-linearly dependent on the added salt concentration. The simulated Setschenow coefficients depend on the DPD bead–bead repulsion amplitudes, and matching to the experimentally determined values provides a systematic method to parameterize the interactions between salt ion beads and surfactant beads. The optimized ion-specific interaction parameters appear to be transferrable and follow the same trends as the empirical Hofmeister series.


■ INTRODUCTION
The presence of ions in a solution affects numerous chemical and biophysical phenomena. The magnitude of these effects often follow ion-specific trends such as the Hofmeister series, 1 which was initially introduced to rank the propensity of salts to decrease the solubility of proteins (salting out) and was subsequently discovered to hold for other phenomena such as partitioning between two liquid phases, 2,3 macromolecular conformational transitions, 4 enzyme activity, 5 protein denaturation, 6,7 viscosity, and critical micelle concentrations (CMCs) of surfactant solutions. 8−10 For a nonelectrolyte, the salting out effect is described by the Setschenow equation 11 f kC ln s salt = − (1) which describes how the activity coefficient f of an uncharged solute depends on the salt concentration C salt . In this, k s is an empirical, salt-specific Setschenow coefficient. We note that the Setschenow relation can be expressed in terms of natural or base 10 logarithms; for this work, we use the natural logarithm "ln". Also, regarding the level of accuracy to which we are working, it is not necessary to distinguish between molar and molal salt concentrations in eq 1 since the difference typically amounts to only a few percent for concentrations less than 1 M at room temperature and pressure. Moreover, since the Setschenow coefficients are arguably defined by their limiting values in eq 1 as C salt tends to zero, this obviates the need to consider the densities of the salt solutions.
While the accurate determination of Setschenow coefficients relies on experimental methods, several models have been developed to predict the effect of a particular salt on a molecule in solution. Among these, an early attempt was made by Debye and MacAulay; 12 then, McDevit and Long later correlated k s for benzene, as a solute, to the change in the volume of the solvent when salt is added; 13 Conway et al. used dielectric saturation to extend the theory to polyions; 14 and Masterton and Lee adopted a scaled particle theory to derive a general expression for the salt effect on benzene derivatives. 15 The use of empirical parameters and the relatively low accuracy of these models limit their use to only a few systems. 16 Wen-Hui et al. linked the Setschenow coefficient for sodium chloride to the Le Bas volume (V LB ) through a simple linear correlation, k s = 0.0018 V LB , 17 while Gould used the intrinsic solubility of the solute. 18 Ni et al. showed how a linear correlation with the partition coefficient of the solute outperformed previous models, 19,20 but data availability limited the study to sodium chloride. More recently, Zhou 21 developed an electrostatic theory to describe the interaction of macromolecules with salt ions, showing good correlation with experimental data for protein solubility and stability. The salting out effect can be directly linked to an increase in surface tension, and Li et al. studied the change in solvation free energy of small molecules in the presence of different salts, concluding that Setschenow coefficients can be explained by the formation of nonpolar cavities in the salt solution and are not due to the direct interaction between solutes and ions. 22 For surfactant solutions, the use of salts to control the CMC, micelle size and shape, and the correlated viscosity makes them an important formulation adjunct for many industrial applications, such as home and personal care products. It is possible to rewrite eq 1 to describe the salt effect on the CMC of nonionic surfactants as 23 where CMC and CMC 0 are, respectively, the CMCs in aqueous salt solution and in pure water. An attractive approach to surfactant formulation is to use computer simulations to supplement or even replace approximate theories and laborious experimental studies. Molecular dynamics (MD) has been extensively used to provide insights into the molecular details, for example, for protein−ion interactions of a single protein in a salt solution. 24−26 Studies have been carried out on the effects of salt on the water structure, 27,28 on the thermodynamics of hydration, 29 and on molecular association. 21,30 For example, Thomas and Elcock used MD to explore the correlation between the hydrophobic effect and water−ion hydrogen bonding using the Hofmeister series. 31 However, statistically meaningful sample sizes are required, and the small size of the systems accessible using MD can be a problem for uncertainty quantification. 32 Due to the high computational costs of MD studies of surfactants, coarse-grained (CG) approaches are more usually deployed to explore micelle formation in surfactant solutions. 33 Salts have been incorporated in CGMD for both hard core and soft core repulsion methods. 34−41 The MARTINI force field describes hard core repulsion using a shifted Lennard-Jones potential, where the parameters depend on the bead type and are optimized using experimental solubility data. 42,43 Alternative methods for parameterization of soft core repulsion have recently been proposed. 37−39 Here, we use dissipative particle dynamics (DPD), a soft core CGMD methodology, which has been developed quite extensively for surfactant simulations. DPD was first introduced by Hoogerbrugge and Koelman 44 and later modified by Espanõl and Warren 45 to satisfy Gibbs− Boltzmann statistics in a canonical NVT ensemble. 46 In DPD, a surfactant solution is modeled using soft particles called DPD beads. The solvent is represented by beads that represent two water molecules, salt ions are represented by adding charges to some of the water beads, and surfactant molecules are represented by a collection of connected DPD beads, which represent the various chemical subgroups. Dissipative and random forces between DPD beads provide a pairwise momentum-conserving thermostat, 45 but the main pairwise interactions are soft, short-range repulsions derived from i k j j j j j j y where a ij is the amplitude of the interaction between beads i and j, r ij is the distance between the two beads, and R ij represents the range of the interaction (cutoff distance). It is through the a ij and R ij parameters that chemical specificity is captured, and recent systematic approaches provide transferable DPD force fields for which the same set of parameters and fragmentation strategy can be used to generate a bead representation of different molecules in different environments.
For this work, we use a recent DPD force field which has been extensively validated for surfactant simulations. 47,48 In the case of added salt, as well as for ionic surfactants, electrostatic interactions must also be included. Methods to implement electrostatics in DPD have received much attention in recent years, with the most common approach being the Ewald method. 33,38,49−53 However, in most studies, salt ions are represented simply as charged water beads, 37,54,55 and very little work has been done to parameterize the (nonelectrostatic) short-range DPD repulsions with these charged beads despite this being an obvious target to capture ionspecific trends such as the Hofmeister series. Along these lines, Mayoral and Nahmad-Achar developed a parameterization for the repulsion amplitudes a ij for the charged beads used in DPD based on the dependency of experimental Flory−Huggins χ parameters on salt concentration. 56 More recently, Nieto-Draghi and Rousseau proposed a parameterization procedure for electrolytes in an aqueous solution based on osmotic pressure, 57 but only ion−ion and ion−water bead interactions were investigated.

■ APPROACH
A systematic strategy to parameterize the interactions between salt beads and the DPD beads which represent the surfactant chemistry has not yet been proposed. Here, we present such a method based on matching to experimental Setschenow coefficients. First, we identify the key interactions between salts and surfactants which affect the CMC. Second, we demonstrate that DPD simulations of surfactants show the same log-linear dependence of the CMC on salt concentration as that reported in eq 2, which allows us to define the DPD equivalents to the empirical Setschenow coefficients. In the third step, we correlate these Setschenow coefficients with the repulsion parameters. Then, in the final step, we utilize the DPD length-scale mapping to match with experimentally determined Setschenow coefficients, providing a systematic basis for fixing the DPD interaction parameters between ions and surfactant beads.
Speciation and micelle formation in aqueous solutions is one of the most common targets of DPD surfactant simulations. 33,58 As already described, treating ions as charged water beads is a simple and convenient approach that has been used to obtain qualitatively good results for several systems. 38,39 However, this approach assumes that all ions are the same, which is in contrast to the experimental evidence for ionspecific trends. Table 1 shows that the experimental Setschenow coefficients that describe the effects of different salts on four nonionic surfactants, n-octanoyl-N-methylglucamine (MEGA8), n-nonanoyl-N-methylglucamine (MEGA9), n-octyl glucopyranoside (GLUCO8), and hexaethylene glycol monododecyl ether (C12E6), cover a wide range of values.
We selected two nonionic surfactants to use for initial parameterization. 60,61 DPD simulations perform most reliably for CMC values in the range 1−100 mM: if the CMC is too high, there is no distinction between submicellar aggregates and true micelles because of the overlap between the submicellar and micellar populations, and conversely, if the CMC is too low, the small number of free surfactants in solution leads to uncertainty. 38,47,62 Accordingly, we selected GLUCO8 and MEGA9, and the CG representations of these surfactants are shown in Figure 1. Table 1 includes Setschenow coefficients for some divalent ions (calcium, monohydrogen phosphate, sulfate, and carbonate). We include these for completeness since they may be useful in future work. To avoid complications arising from the nonideality of salt The Journal of Physical Chemistry B pubs.acs.org/JPCB Article solutions containing multivalent ions, which may be significant, the present study focuses on the monovalent ions in Table 1.

■ SIMULATION DETAILS
DPD simulations were run for all the surfactant systems at different salt concentrations. The CG representation of the surfactants and the repulsion parameters for the beads were obtained from previous studies ( Figure 1 and Table 2). 47,48 The cutoff distances in Table 2 were assigned as given in Anderson et al. 63 The repulsion amplitude between water beads is chosen conventionally as 25 k B T so that the pressure of pure water in DPD units is 23.7. 63 In the model, the cutoff distance R ij between water beads is defined as the DPD unit of length r c , pure water is represented by water beads at a reduced density ρr c 3 = 3, and we suppose that each water bead represents N m water molecules, where N m is the so-called (water bead) mapping number. 64 To fit with the chosen DPD force field, we use N m = 2. 45,62 If V m ≈ 18 × 10 −6 m 3 is the molar volume of water, one can deduce that the volume of 1 mol of DPD volume elements N A r c 3 = ρr c 3 × N m × V m ≈ 0.108 × 10 −3 m 3 ≈ 0.108 litres, and hence, r c = 5.64 Å. The DPD volume can be used to convert the number of salt beads (N salt ) to molar concentration units; namely, if the simulation box side is L, then Molecularly bonded beads are held together with a harmonic spring potential where k b = 150 k B T, and a three-body angular potential where k a = 5 k B T. The equilibrium distance r 0 and the equilibrium angle θ 0 were assigned using the method previously reported. 48 Simulations were performed in a cubic box of side L = 30 r c with the total number of beads equal to 81,000. They were run   for 4 × 10 6 steps with a time step of 0.01 in DPD time units, starting from a random initial configuration. By measuring the diffusion of small molecules in a related DPD model with a comparable level of coarse graining, Sevink and Fraaije determined the underlying DPD time unit ≈ 50 ps, so the time step in our simulations should correspond to about 0.5 ps, 65 making the total simulation run time 1−2 μs. This timescale is long enough for micelles to form and to equilibrate. Simulations were run in the isothermal−isobaric ensemble (NPT) using the standard velocity Verlet integration algorithm. 66 Trajectory files were collected every 10 3 time steps. Simulations were performed using DL_MESO (version 2.7). 67 Post-simulation trajectory analysis was performed using a combination of the UMMAP analysis tool 68 and bespoke analysis scripts. DPD simulations were run at 4, 5, and 6 wt % for all surfactants. The aggregation number distribution is a plot of population P(N) versus aggregation number (N), and this distribution can be used to discriminate between monomers and submicellar aggregates (designated free surfactant) and stable micelles as described in previous studies. 47,48 By plotting this distribution, one discerns a region depleted in stable micelles, which allows the definition of a value of N cut to separate the free surfactants (N < N cut ) from micelles (N > N cut ). For each simulation, the minimum in P(N) from the aggregation number distribution was used as N cut . The CMC was calculated as the total concentration of free surfactants after reaching the equilibrium (typically after 5 × 10 5 steps). The values of CMC did not vary significantly with surfactant concentration, and the average values are quoted.
As a starting point, each salt ion was represented by a positive bead or a negative bead, with the same R ij and a ij as that of water (i.e., as charged water beads), and then, the a ij values were varied to study the effect on the calculated CMC of the surfactant. The Slater-type charge smearing for the electrostatic interactions proposed by Gonzaĺez-Melchor et al. was adopted. 69 For a pair of particles, the electrostatic interaction potential can be written as where r ij is the distance between particles i and j, q i and q j are the ion valences, and β* is the Slater smearing parameter (set to be equal to 0.929 r c

■ RESULTS AND DISCUSSION
First, the effect of changing the values of a ij for interactions between salt ions and other beads was investigated though the effect on the CMC. Figure 2 shows the relationship between the CMC and the a ij parameters for ion−water (blue), ion−tail (green), and ion−head group (yellow) interactions. Beads T, T2, and C2 were considered part of the hydrophobic tail, and beads AM3, OH, OH′, OH2, and AC′ were considered part of the hydrophilic head group. Figure 2 compares the results for MEGA9 and GLUCO8 obtained using a salt concentration of 1 M with the experimental values (dotted line). The plots show that the ion−tail interactions have a much bigger effect in decreasing the CMC values than the ion−head and ion−water interactions and that the best way to approach values closer to the experimental values is by further increasing the repulsion parameter for the ion−tail interactions. Therefore, we focused our attention on these interaction parameters, which is in line with the conclusions obtained by Mukerjee in his studies on C x E y surfactants, where the contribution of the hydrophilic head to k s for neutral surfactants was reported to be negligible. 23 Ion−tail a ij values were screened from 35 k B T to 150 k B T for both the anion (a ij − ) and cation (a ij + ), and the CMC values were calculated for both surfactants. The result for GLUCO8 with a ij + = a ij − reported in Figure 3a shows a nonlinear relationship, which reaches a plateau around a ij = 140 k B T. The CMC value increases sharply for small a ij values, which is consistent with the negative k s values observed for polyoxyethylenes. 71 The effect of salt concentration on the CMC value was then studied. The CMC values were averaged from three independent simulation runs with surfactant concentrations of 4, 5, and 6 wt %. These simulations were repeated using different ion−tail a ij values and for a range of salt   Figure 3b shows an example of the results. In all cases, a clear log-linear relationship was found between CMC and salt concentration. Slopes of plots of the CMC values obtained in DPD simulations versus salt concentration were used to calculate values of k s as a function of the ion−tail repulsion parameters used in the simulations. The relationship between the calculated values of k s and the values of a ij − and a ij + for the anion−tail and cation−tail interactions is shown in Figure 4 for GLUCO8 and MEGA9. The surfaces colored in green in Figure 4 are the best fit for the following relationship i k j j j j j j y where A and B are constants that depend on the surfactant ( Table 3).
The interesting property of eq 6 is that the terms that describe interactions with the anion and the cation appear separately so that a salt can be described simply as the sum of the individual effects of the anion and the cation on the CMC.
Having established values of A and B for GLUCO8 and MEGA9, the experimentally determined values of k s measured for these surfactants in different salt solutions can be used in conjunction with eq 6 to derive the repulsion parameters required to describe the individual ions. The ion−tail repulsion parameters for lithium, sodium, and potassium cations and for chloride, bromide, and iodide anions were obtained by fitting the experimental values of k s for GLUCO8 and MEGA9 in Table 1 to eq 6 using a generalized reduced gradient nonlinear method. Having established repulsion parameters for a range of different anions and cations, the experimentally determined values of k s measured for different surfactants in the corresponding salt solutions can be used in conjunction with eq 6 to derive the constants A and B required to describe the surfactants. The parameters for MEGA8 and C12E6 were obtained by fitting the experimental values of k s in Table 1 to eq 6, and the results are shown in Table 3. The repulsion parameters for the ions listed in Table 4 were then optimized by using the experimental values of k s for all four surfactants in eq 6 and the values of A and B in Table 3. The parameters in Tables 3 and 4 provide an excellent description of the experimental k s values for all four surfactants in solutions of 10 different salts ( Figure 5). The results suggest that eq 6 can be used to predict the effects of salts on other surfactants, provided sufficient experimental data are available to estimate the A and B parameters for the surfactant ( Table 5).
The parameters for thiocyanate and caesium ions in Table 4 were obtained from single measurements and should be considered less reliable than the other values. The parameter for the fluoride ion could not be determined by changing the ion−tail repulsion parameter because the a ij value required was too large. It is clear that fluoride has the highest repulsion parameter of all of the anions, but some other interactions must be involved to account for the behavior of these systems, such as a specific interaction with the surfactant or a direct effect on the counterion. The trends in a ij values are F − > Cl − > Br − > NO 3 − > I − > SCN − for anions and Na + > K + > Cs + > Li + for cations. With the exception of K + , which precedes Na + in     Our observation that the Setschenow trends can be captured in the present model by tuning the ion−tail interactions supports the notion that the Setschenow coefficients are a measure of the solvent "quality" in these systems. Thus, our results underscore the idea that, in this context, the Hofmeister series reflects changes in the hydrophobic effect acting on the surfactant tails. 73 We emphasize in this respect that apart from the ion−tail interactions in Table 4, the ions are otherwise treated as charged water beads (cf., the "W" beads given in Table 2). This means for instance that trends in ion activities are not reproduced in the present model. This could be solved by combining the approach with, for example, the parameterization strategy for ion−ion and ion−water interactions proposed by Nieto-Draghi and Rousseau, 57 which could then be extended to include multivalent ions. Further work will also be required to extend the approach to ionic surfactants and to test the effect of changes in the nature of the hydrocarbon tail. The case of polyions (polyelectrolytes) is also of considerable interest; however, these are often involved with very specific effects such as adsorption onto surfactant micelles 74 and would require separate treatment.